#include "sourceField.hpp"
#include "fieldFunctions.hpp"
#include "util.hpp"
#include "typeHelpers.hpp"
#include "field.hpp"
#include <cassert>

namespace hms
{

#define SFHEAD(RET) \
template<Uniformity isUniform, typename MeshType> \
RET SourceField<isUniform, MeshType>::

SFHEAD()
SourceField(
	const MeshBase<MeshType>& mesh,
	ArrayXXs&& vals
) :
	FieldBase<Type>(mesh)
{
	values( std::move(vals) );
	m_useInterpolation = false;
}

SFHEAD()
SourceField(
	const MeshBase<MeshType>& mesh,
	std::vector<std::pair<scalar,ArrayXXs>>&& vals,
	const Time& t,
	bool useInterpolationArg
) :
	FieldBase<Type>(mesh)
{
	timeSeries( std::move(vals) );
	time(t);
	m_useInterpolation = useInterpolationArg;
}

SFHEAD()
SourceField(
	const MeshBase<MeshType>& mesh,
	std::vector<FunctionParserWrapper>&& parsers,
	const Time& t,
	const Field<MeshType>& z
) :
	FieldBase<Type>(mesh)
{
	functions( std::move(parsers) );
	time(t);
	zField(z);
}

SFHEAD()
SourceField(
	const MeshBase<MeshType>& mesh,
	std::function<void(ArrayXXs&)>&& func,
	Index nComponents
) :
	FieldBase<Type>(mesh)
{
	callable( std::move(func), nComponents );
}

SFHEAD(void)
values( ArrayXXs&& vals ){
	assertFieldValuesSize<isUniform>( vals, this->mesh() );
	m_values = std::move(vals);
	if constexpr (!isUniform)
		setCorners(m_values, this->mesh(), FieldBase<Type>::cornerValue);
	m_updateMethod = UpdateMethod::none;
	m_useInterpolation  = false;
}

SFHEAD(void)
timeSeries( std::vector<std::pair<scalar,ArrayXXs>>&& vals ){
	#ifndef NDEBUG
	for ( const auto& pair : vals ){
		assertFieldValuesSize<isUniform>( pair.second, this->mesh() );
	}
	assert( std::is_sorted(vals.cbegin(), vals.cend(),
		[](const auto& pair1, const auto& pair2){
			return pair1.first < pair2.first;
		} )
	);
	#endif
	if ( vals.size() > 1 ){
		m_updateMethod = UpdateMethod::timeSeries;
		m_timeSeries = std::move(vals);
		m_values = m_timeSeries[0].second;
		if constexpr (!isUniform){
			for ( auto& pair : m_timeSeries )
				setCorners(pair.second, this->mesh(), FieldBase<Type>::cornerValue);
		}
	} else
		values( std::move( vals[0].second ) );
	
}

SFHEAD(void)
functions( std::vector<FunctionParserWrapper>&& parsers ){
	m_functions = std::move(parsers);
	m_updateMethod = UpdateMethod::function;
	m_useInterpolation  = false;
	resizeFieldValues<isUniform>( m_values, m_functions.size(), this->mesh() );
	if constexpr (!isUniform)
		setCorners(m_values, this->mesh(), FieldBase<Type>::cornerValue);
}

SFHEAD(void)
callable( std::function<void(ArrayXXs&)>&& func, Index nComponents ){
	m_updater = std::move(func);
	m_updateMethod = UpdateMethod::callable;
	m_useInterpolation  = false;
	resizeFieldValues<isUniform>( m_values, nComponents, this->mesh() );
}

SFHEAD(void)
time( const Time& t ){

	m_time = std::experimental::make_observer( &t );

	if ( m_updateMethod == UpdateMethod::timeSeries ){
		assert( !m_timeSeries.empty() && "Initialise time series first!");
		auto iter { findValuesAtTime(m_timeSeries, time().current()) };
		if ( iter == m_timeSeries.cend() )
			m_currentIndex = m_timeSeries.size() - 1;
		else
			m_currentIndex = iter - m_timeSeries.cbegin();
	}
}


SFHEAD(void)
zField( const Field<MeshType>& z ){
	if constexpr (!isUniform){
		assert( !z.empty() );
		m_zField = std::experimental::make_observer( &z );
	}
}



SFHEAD(auto)
storage() const -> const ArrayXXs& { return m_values; }

SFHEAD(auto)
storage() -> ArrayXXs& { return m_values; }


SFHEAD(auto)
timeSeries() const -> const std::vector<std::pair<scalar,ArrayXXs>>& {
	assert( m_updateMethod == UpdateMethod::timeSeries );
	return m_timeSeries;
}

SFHEAD(const std::vector<FunctionParserWrapper>&)
functions() const {
	assert( m_updateMethod == UpdateMethod::function );
	return m_functions;
}

SFHEAD(const Time&)
time() const {
	assert(m_time && "Source field is missing time object!");
	return *m_time;
}


SFHEAD(auto)
updateMethod() const -> UpdateMethod {
	return m_updateMethod;
}


SFHEAD(bool)
useInterpolation() const {
	return m_useInterpolation;
}

SFHEAD(void)
useInterpolation(bool val) {
	assert( m_updateMethod == UpdateMethod::timeSeries );
	m_useInterpolation = val;
}

SFHEAD(void)
updateValues(){
	using U = UpdateMethod;
	switch (m_updateMethod){
		case U::none:
			break;
		case U::timeSeries:
			updateValues_timeSeries_preventNesting();
			break;
		case U::function:
			updateValues_function();
			break;
		case U::callable:
			updateValues_callable();
			break;
	}
}

SFHEAD(void)
updateValues_timeSeries(){
	/* increment current index if necessary */
	bool indexChanged {false};
	using SizeType = remove_qualifiers<decltype( m_timeSeries.size() )>;
	while (
		static_cast<SizeType>(m_currentIndex+1) < m_timeSeries.size() &&
		time().current() >= m_timeSeries[m_currentIndex+1].first
	){
		++m_currentIndex;
		indexChanged = true;
	}
	/* current time values */
	const auto& currentPair { m_timeSeries[m_currentIndex] };

	/* if the index is pointing to the last element already, no later time
	 * values are available */
	bool atEnd {
		static_cast<SizeType>(m_currentIndex+1) == m_timeSeries.size()
	};

	if ( m_useInterpolation && !atEnd ){
		const auto& nextPair { m_timeSeries[m_currentIndex+1] };
		scalar ratio {
			(time().current() - currentPair.first) /
			( nextPair.first - currentPair.first)
		};
		m_values = hms::interpolate(
			currentPair.second, nextPair.second, ratio
		);
	} else if (indexChanged) {
		m_values = currentPair.second;
	}
}

SFHEAD(void)
updateValues_timeSeries_preventNesting(){
	if ( omp_in_parallel() ){
		OMP_PRAGMA(single nowait)
			updateValues_timeSeries();
	} else {
		updateValues_timeSeries();
	}
}

SFHEAD(void)
updateValues_function(){
	/* is parallelised */
	setFieldByFunctions<isUniform>(
		m_values, m_functions, this->mesh(), time(),
		m_zField
	);
}

SFHEAD(void)
updateValues_callable(){
	/* updater must be parallelised! */
	m_updater(m_values);
	assertFieldValuesSize<isUniform>( m_values, this->mesh() );
}


SFHEAD(auto)
valuesImp() const -> ConstValueType {
	assert( m_values.rows() != 0 && m_values.cols() != 0 );

	if constexpr (isUniform){
		return m_values.replicate( 1, nFieldCols( this->mesh() ) );
	} else {
		return m_values;
	}
}

SFHEAD(auto)
zField() const -> const Field<MeshType>& {
	assert( m_zField && "zField pointer is nullptr!" );
	return *m_zField;
}

#undef SFHEAD

#define INSTANTIATE_FOR_ALL_MESHTYPES(MeshType) \
template class SourceField<Uniform, MeshType>; \
template class SourceField<NonUniform, MeshType>;
#include "meshInstantiator.hpp"

} // namespace hms
